Implementation of the reciprocal square root in MPFR

We describe the implementation of the reciprocal square root — also called inverse square root — as a native function in the MPFR library. The difficulty is to implement Newton’s iteration for the reciprocal square root on top’s of GNU MP’s mpn layer, while guaranteeing a rigorous 1/2 ulp bound on the roundoff error. The reciprocal square root is an important function in 3D graphics, for the normalization of 3D vectors, and as such has received much attention in the literature [5]. Indeed, given (x, y, z), compute t = x2 + y2 + z2, u = t−1/2, x′ = ux, y′ = uy, z′ = uz, then (x′, y′, z′) is a normalized vector with x′2 + y′2 + z′2 = 1. A nice implementation attributed to Greg Walsh by Wikipedia is the following: float invSqrt(float x) { float xhalf = 0.5f*x; int i = *(int *)&x; i = 0x5f3759df (i >> 1); x = *(float *)&i; x = x * (1.5f xhalf * x * x); return x; } See also IBM’s US patent 6654777 (Single precision inverse square root generator) issued on November 25, 2003. The reciprocal square root is an interesting function per se, since it can be computed with Newton’s iteration without any division, similarly to the reciprocal. Karp and Markstein have shown how the square root can be efficiently computed, by modifying the last Newton iteration of the reciprocal square root [4]. In the context of the MPFR library [2], our aim is to provide a native implementation of the reciprocal square root which guarantees correct rounding. By “native implementation”, we mean it should rely on the GMP mpn layer only, which provides operations on arrays of words (addition, subtraction, multiplication) [3]. MPFR reciprocal square root (function mpfr rec sqrt) is based on Ziv’s strategy and the mpfr mpn rec sqrt function, which given a precision p, and an input 1 ≤ A < 4, returns an approximation x satisfying x− 1 2 · 2−p ≤ A−1/2 ≤ x+ 2−p. The mpfr mpn rec sqrt function is based on Newton’s iteration and the following lemma [1]: Date: March 2008. Centre de Recherche INRIA Nancy Grand Est, Projet CACAO Bâtiment A, 615 rue du Jardin Botanique, F-54602 Villers-lès-Nancy Cedex, [email protected].